home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_5 / a5_2.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.4 KB  |  135 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 5.2 (Least Squares Polynomial).
  9. % Section    5.2,    Curve Fitting, Page 278
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % This program finds the least squares polynomial Pm(x), given
  13.  
  14. % a set of data points { (x , y ), (x , y ) ,..., (x , y ) }.
  15. %                          1   1     2   2          n   n
  16.  
  17. % The abscissas and ordinates are stored in X and Y, respectively.
  18.  
  19. % X = [x , x  ,..., x ];  Y = [y , y  ,..., y ];
  20. %       1   2        n          1   2        n
  21.  
  22. % Remark. lspoly.m is used for Algorithm 5.2
  23.  
  24. pause % Press any key to continue.
  25.  
  26. clc;
  27. %    Place the abscissas for the points in  X.
  28.  
  29. %    Place the ordinates for the points in  Y.
  30.  
  31. % Place the degree of the polynomial in  m.
  32.  
  33. X = [-3  -2  -1   0  1   2   3  4];
  34.  
  35. Y = [ 3   2  1.5  1  1  1.5  3  5];
  36.  
  37. m = 2;
  38.  
  39. C = lspoly(X,Y,m);
  40.  
  41. pause    % Press any key to find the least squares polynomial.
  42.  
  43. clc; clg;
  44. a = -4;
  45. b = 5;
  46. c = 0;
  47. d = 6;
  48. h = (b-a)/150;
  49. Xs = a:h:b;
  50. Ys = polyval(C,Xs);
  51. axis([a b c d]);...
  52. plot(X,Y,'or',Xs,Ys,'-g');...
  53. hold on;...
  54. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  55. xlabel('x');...
  56. xlabel('x');...
  57. ylabel('y');...
  58. Mx1 = 'Least squares polynomial: y = P';...
  59. Mx2 = [Mx1,num2str(m),'(x).'];...
  60. title(Mx2);...
  61. grid;...
  62. axis;...
  63. hold off;...
  64. shg; pause    % Press any key to continue.
  65.  
  66. points1 = [X;Y]; format short;
  67. Mx1='y = P(x) = c(1)x^m + c(2)x^m-1 +...+ c(m)x + c(m+1)';
  68. Mx2=['A polynomial of degree m = ',num2str(m),' has been fit.'];
  69. Mx3='The coefficients are stored in the array  C = ';
  70. clc,echo off,diary output,...
  71. disp(''),disp(Mx1),disp(Mx2),disp(Mx3),...
  72. disp(''),disp(C'),disp('The given x-y points:'),...
  73. disp('      x         y'),disp(points1'),diary off,echo on
  74. pause    % Press any key to analyze the results.
  75.  
  76. points2 = [X;Y;polyval(C,X);Y-polyval(C,X)]';
  77. Mx4='    x(k)      y(k)      P(x(k))   error';
  78. clc,echo off,diary output,...
  79. disp(''),disp(Mx4),disp(points2),diary off,echo on
  80. pause    % Press any key to continue.
  81.  
  82. clc;
  83. %    Place the abscissas for the points in  X.
  84.  
  85. %    Place the ordinates for the points in  Y.
  86.  
  87. % Place the degree of the polynomial in  m.
  88.  
  89. X = [-3  -2  -1   0  1   2   3  4];
  90.  
  91. Y = [ 3   2  1.5  1  1  1.5  3  5];
  92. m = 3;
  93.  
  94. C = lspoly(X,Y,m);
  95.  
  96. pause    % Press any key to find the least squares polynomial.
  97.  
  98. clc; clg;
  99. a = -4;
  100. b = 5;
  101. c = 0;
  102. d = 6;
  103. h = (b-a)/150;
  104. Xs = a:h:b;
  105. Ys = polyval(C,Xs);
  106. axis([a b c d]);...
  107. plot(X,Y,'or',Xs,Ys,'-g');...
  108. hold on;...
  109. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  110. xlabel('x');...
  111. xlabel('x');...
  112. ylabel('y');...
  113. Mx1 = 'Least squares polynomial: y = P';...
  114. Mx2 = [Mx1,num2str(m),'(x).'];...
  115. title(Mx2);...
  116. grid;...
  117. axis;...
  118. hold off;...
  119. shg; pause    % Press any key to continue.
  120.  
  121. points1 = [X;Y]; format short;
  122. Mx1='y = P(x) = c(1)x^m + c(2)x^m-1 +...+ c(m)x + c(m+1)';
  123. Mx2=['A polynomial of degree m = ',num2str(m),' has been fit.'];
  124. Mx3='The coefficients are stored in the array  C = ';
  125. clc,echo off,diary output,...
  126. disp(''),disp(Mx1),disp(Mx2),disp(Mx3),...
  127. disp(''),disp(C'),disp('The given x-y points:'),...
  128. disp('      x         y'),disp(points1'),diary off,echo on
  129. pause    % Press any key to analyze the results.
  130.  
  131. points2 = [X;Y;polyval(C,X);Y-polyval(C,X)]';
  132. Mx4='    x(k)      y(k)      P(x(k))   error';
  133. clc,echo off,diary output,...
  134. disp(''),disp(Mx4),disp(points2),diary off,echo on
  135.